#<(GPL, copyleft, https://www.gnu.org/licenses/gpl-3.0.de.html)>
<This code shows scalewise confirmatory factor analyses.>
Copyright (C) <2021>  <Walter, B.; Zehtner, R. I. & Hermann, A.>
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

#This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

#You should have received a copy of the GNU General Public License
along with this program.  If not, see 
https://www.gnu.org/licenses/ <https://www.gnu.org/licenses/>

#================================================================

#analyses were computed with R version 4.0.4 (R Core Team, 2019) 

#To run this program you need to specify the path of data as well as the data frames
#Note: Items 19 and 28 were deleted before running the CFAs

#================================================================

rm(list = ls())
setwd('D:/FEQ_GR/Studie1/Berechnungen')


library(tidyverse)
library(foreign)
library(psych)
library(psychTools)
library(lattice)

library(lavaan) # for CFA
library(semTools) # for additional functions in SEM
library(semPlot) # for path diagram

library(MBESS) # for CIs of reliability estimates

library(RGenData)



# ++++++++++++++++++++++
# create data frames
# ++++++++++++++++++++++
df1 <- read.spss( "data_study1_Subsample1.sav", use.value.labels=FALSE, to.data.frame=TRUE )
names(df1)
feq_sample_1 <- df1[36:73]

df2 <- read.spss( "data_study1_Subsample2.sav", use.value.labels=FALSE, to.data.frame=TRUE )
names(df2)
feq_sample_2 <- df2[36:73]

feq_total <- bind_rows(feq_sample_1, feq_sample_2)

str(feq_total)

#############################
#CFAs with 2 and 3 factors as used in literature

# f2 scale: Positive expressiveness (P)

# unidimensionality
fa2_model_scaleP <- "
fa2_scaleP =~ FEQ_01 + FEQ_06 + FEQ_16 + FEQ_17 + FEQ_18 + FEQ_23 + FEQ_26 + FEQ_33 + 
              FEQ_39 + FEQ_02 + FEQ_03 + FEQ_13 + FEQ_21 + FEQ_22 + FEQ_30 + FEQ_31 + 
              FEQ_35 + FEQ_38 + FEQ_40
"

s1.fa2.cfa_scaleP <- cfa(fa2_model_scaleP, data = feq_sample_1, estimator = "MLM", mimic = "Mplus")

summary(s1.fa2.cfa_scaleP, fit.measures = T, standardized = T)


# reliability with Item FEQ_13
s1.fa2_scaleP <- feq_sample_1 %>% select(FEQ_01, FEQ_06, FEQ_16, FEQ_17, FEQ_18, FEQ_23, FEQ_26, FEQ_33, FEQ_39, 
                                         FEQ_02, FEQ_03, FEQ_13, FEQ_21, FEQ_22, FEQ_30, FEQ_31, FEQ_35, FEQ_38, FEQ_40)

set.seed(42)
ci.reliability(data = s1.fa2_scaleP, type="hierarchical", conf.level = 0.95, interval.type="bca", B=1000)
alpha.fa2_scaleP = alpha(s1.fa2_scaleP)
print(alpha.fa2_scaleP, digits = 3)

################################delete item 13

# unidimensionality without item 13
fa2_model_scaleP <- "
fa2_scaleP =~ FEQ_01 + FEQ_06 + FEQ_16 + FEQ_17 + FEQ_18 + FEQ_23 + FEQ_26 + FEQ_33 + 
              FEQ_39 + FEQ_02 + FEQ_03 + FEQ_21 + FEQ_22 + FEQ_30 + FEQ_31 + 
              FEQ_35 + FEQ_38 + FEQ_40
"

s1.fa2.cfa_scaleP <- cfa(fa2_model_scaleP, data = feq_sample_1, estimator = "MLM", mimic = "Mplus")

summary(s1.fa2.cfa_scaleP, fit.measures = T, standardized = T)


# reliability without item 13
s1.fa2_scaleP <- feq_sample_1 %>% select(FEQ_01, FEQ_06, FEQ_16, FEQ_17, FEQ_18, FEQ_23, FEQ_26, FEQ_33, FEQ_39, 
                                         FEQ_02, FEQ_03, FEQ_21, FEQ_22, FEQ_30, FEQ_31, FEQ_35, FEQ_38, FEQ_40)

set.seed(42)
ci.reliability(data = s1.fa2_scaleP, type="hierarchical", conf.level = 0.95,
               interval.type="bca", B=1000)
alpha.fa2_scaleP = alpha(s1.fa2_scaleP)
print(alpha.fa2_scaleP, digits = 3)

##################################Crossvalidation in subsample 2 (without item 13)

# unidimensionality
fa2_model_scaleP <- "
fa2_scaleP =~ FEQ_01 + FEQ_06 + FEQ_16 + FEQ_17 + FEQ_18 + FEQ_23 + FEQ_26 + FEQ_33 + 
              FEQ_39 + FEQ_02 + FEQ_03 + FEQ_21 + FEQ_22 + FEQ_30 + FEQ_31 + 
              FEQ_35 + FEQ_38 + FEQ_40
"

s2.fa2.cfa_scaleP <- cfa(fa2_model_scaleP, data = feq_sample_2, estimator = "MLM", mimic = "Mplus")

summary(s2.fa2.cfa_scaleP, fit.measures = T, standardized = T)


# reliability ohne Item FEQ_13
s2.fa2_scaleP <- feq_sample_2 %>% select(FEQ_01, FEQ_06, FEQ_16, FEQ_17, FEQ_18, FEQ_23, FEQ_26, FEQ_33, FEQ_39, 
                                         FEQ_02, FEQ_03, FEQ_21, FEQ_22, FEQ_30, FEQ_31, FEQ_35, FEQ_38, FEQ_40)

set.seed(42)
ci.reliability(data = s2.fa2_scaleP, type="hierarchical", conf.level = 0.95,
               interval.type="bca", B=1000)
alpha.fa2_scaleP = alpha(s2.fa2_scaleP)
print(alpha.fa2_scaleP, digits = 3)

#############################################
## f2 scale: Negative expressiveness (N)

# unidimensionality
fa2_model_scaleN <- "
fa2_scaleN =~ FEQ_04 + FEQ_05 + FEQ_07 + FEQ_09 + FEQ_11 + FEQ_12 + FEQ_24 + FEQ_27 + 
              FEQ_36 + FEQ_37 + FEQ_08 + FEQ_10 + FEQ_14 + FEQ_15 + FEQ_20 + FEQ_25 + 
              FEQ_29 + FEQ_32 + FEQ_34
"

s1.fa2.cfa_scaleN <- cfa(fa2_model_scaleN, data = feq_sample_1, estimator = "MLM", mimic = "Mplus")

summary(s1.fa2.cfa_scaleN, fit.measures = T, standardized = T)


# reliability
s1.fa2_scaleN <- feq_sample_1 %>% select(FEQ_04, FEQ_05, FEQ_07, FEQ_09, FEQ_11, FEQ_12, 
                                         FEQ_24, FEQ_27, FEQ_36, FEQ_37, FEQ_08, FEQ_10,
                                         FEQ_14, FEQ_15, FEQ_20, FEQ_25, FEQ_29, FEQ_32, FEQ_34)

set.seed(42)
ci.reliability(data = s1.fa2_scaleN, type="hierarchical", conf.level = 0.95,
               interval.type="bca", B=1000)
alpha.fa2_scaleN = alpha(s1.fa2_scaleN)
print(alpha.fa2_scaleN, digits = 3)

#2 items with negative loadings, i.e., 
#++Item 29 (Für zu spät kommen entschuldigen)
#++Item 34 (Sagen, wie verletzt man ist)
#4 items with insuffiecent loadings, i.e., 
#++Item 25, Item 14, Item 20, Item 37 

###########################CFA -6 items Subsample 1
## f2 scale: Negative expressiveness (N)

# unidimensionality minus 6 Items
fa2_model_scaleN <- "
fa2_scaleN =~ FEQ_04 + FEQ_05 + FEQ_07 + FEQ_09 + FEQ_11 + FEQ_12 + FEQ_24 + FEQ_27 + 
              FEQ_36 + FEQ_08 + FEQ_10 + FEQ_15 + FEQ_32
"

s1.fa2.cfa_scaleN <- cfa(fa2_model_scaleN, data = feq_sample_1, estimator = "MLM", mimic = "Mplus")

summary(s1.fa2.cfa_scaleN, fit.measures = T, standardized = T)


# reliability minus 6 Items 
s1.fa2_scaleN <- feq_sample_1 %>% select(FEQ_04, FEQ_05, FEQ_07, FEQ_09, FEQ_11, FEQ_12, 
                                        FEQ_24, FEQ_27, FEQ_36, FEQ_08, FEQ_10, FEQ_15, FEQ_32)

set.seed(42)
ci.reliability(data = s1.fa2_scaleN, type="hierarchical", conf.level = 0.95,
               interval.type="bca", B=1000)
alpha.fa2_scaleN = alpha(s1.fa2_scaleN)
print(alpha.fa2_scaleN, digits = 3)

#############################CFA -6 items in Subsample 2
## f2 scale: Negative expressiveness (N)

# unidimensionality minus 6 Items, Subsample 2
fa2_model_scaleN <- "
fa2_scaleN =~ FEQ_04 + FEQ_05 + FEQ_07 + FEQ_09 + FEQ_11 + FEQ_12 + FEQ_24 + FEQ_27 + 
              FEQ_36 + FEQ_08 + FEQ_10 + FEQ_15 + FEQ_32
"

s2.fa2.cfa_scaleN <- cfa(fa2_model_scaleN, data = feq_sample_2, estimator = "MLM", mimic = "Mplus")

summary(s2.fa2.cfa_scaleN, fit.measures = T, standardized = T)


# reliability minus 6 Items, Subsample 2 
s2.fa2_scaleN <- feq_sample_2 %>% select(FEQ_04, FEQ_05, FEQ_07, FEQ_09, FEQ_11, FEQ_12, 
                                         FEQ_24, FEQ_27, FEQ_36, FEQ_08, FEQ_10, FEQ_15, FEQ_32)

set.seed(42)
ci.reliability(data = s2.fa2_scaleN, type="hierarchical", conf.level = 0.95,
               interval.type="bca", B=1000)
alpha.fa2_scaleN = alpha(s2.fa2_scaleN)
print(alpha.fa2_scaleN, digits = 3)

#######################################################################
####CFA with three scales: P, NS, ND
#Sample1

#P = positive scale remains the same

#Negative-submissive

# unidimensionality scale NS
fa3_model_scaleNS <- "
fa3_scaleNS =~ FEQ_08 + FEQ_10 + FEQ_14 + FEQ_15 + FEQ_20 + FEQ_25 + 
               FEQ_29 + FEQ_32 + FEQ_34
"

s1.fa3.cfa_scaleNS <- cfa(fa3_model_scaleNS, data = feq_sample_1, estimator = "MLM", mimic = "Mplus")

summary(s1.fa3.cfa_scaleNS, fit.measures = T, standardized = T)


# reliability scale NS
s1.fa3_scaleNS <- feq_sample_1 %>% select(FEQ_08, FEQ_10,
                                         FEQ_14, FEQ_15, FEQ_20, FEQ_25, FEQ_29, FEQ_32, FEQ_34)

set.seed(42)
ci.reliability(data = s1.fa3_scaleNS, type="hierarchical", conf.level = 0.95,
               interval.type="bca", B=1000)
alpha.fa3_scaleNS = alpha(s1.fa3_scaleNS)
print(alpha.fa3_scaleNS, digits = 3)

##################################2 Items deleted: Item 29, Item 34
# unidimensionality without Item 29, 34
fa3_model_scaleNS <- "
fa3_scaleNS =~ FEQ_08 + FEQ_10 + FEQ_14 + FEQ_15 + FEQ_20 + FEQ_25 + 
               FEQ_32
"

s1.fa3.cfa_scaleNS <- cfa(fa3_model_scaleNS, data = feq_sample_1, estimator = "MLM", mimic = "Mplus")

summary(s1.fa3.cfa_scaleNS, fit.measures = T, standardized = T)


# reliability scale NS without Item 29, 34
s1.fa3_scaleNS <- feq_sample_1 %>% select(FEQ_08, FEQ_10,
                                          FEQ_14, FEQ_15, FEQ_20, FEQ_25, FEQ_32)

set.seed(42)
ci.reliability(data = s1.fa3_scaleNS, type="hierarchical", conf.level = 0.95,
               interval.type="bca", B=1000)
alpha.fa3_scaleNS = alpha(s1.fa3_scaleNS)
print(alpha.fa3_scaleNS, digits = 3)

###########################Überprüfung NS an Subsample 2
# unidimensionality without Item 29, 34
fa3_model_scaleNS <- "
fa3_scaleNS =~ FEQ_08 + FEQ_10 + FEQ_14 + FEQ_15 + FEQ_20 + FEQ_25 + 
               FEQ_32
"

s2.fa3.cfa_scaleNS <- cfa(fa3_model_scaleNS, data = feq_sample_2, estimator = "MLM", mimic = "Mplus")

summary(s2.fa3.cfa_scaleNS, fit.measures = T, standardized = T)


# reliability scale NS
s2.fa3_scaleNS <- feq_sample_2 %>% select(FEQ_08, FEQ_10,
                                          FEQ_14, FEQ_15, FEQ_20, FEQ_25, FEQ_32)

set.seed(42)
ci.reliability(data = s2.fa3_scaleNS, type="hierarchical", conf.level = 0.95,
               interval.type="bca", B=1000)
alpha.fa3_scaleNS = alpha(s2.fa3_scaleNS)
print(alpha.fa3_scaleNS, digits = 3)

###########################Negative-dominant
# unidimensionality ND
fa3_model_scaleND <- "
fa3_scaleND =~ FEQ_04 + FEQ_05 + FEQ_07 + FEQ_09 + FEQ_11 + FEQ_12 + FEQ_24 + FEQ_27 + 
              FEQ_36 + FEQ_37
"

s1.fa3.cfa_scaleND <- cfa(fa3_model_scaleND, data = feq_sample_1, estimator = "MLM", mimic = "Mplus")

summary(s1.fa3.cfa_scaleND, fit.measures = T, standardized = T)


# reliability ND
s1.fa3_scaleND <- feq_sample_1 %>% select(FEQ_04, FEQ_05, FEQ_07, FEQ_09, FEQ_11, FEQ_12, 
                                         FEQ_24, FEQ_27, FEQ_36, FEQ_37)

set.seed(42)
ci.reliability(data = s1.fa3_scaleND, type="hierarchical", conf.level = 0.95,
               interval.type="bca", B=1000)
alpha.fa3_scaleND = alpha(s1.fa3_scaleND)
print(alpha.fa3_scaleND, digits = 3)

########################ND in Subsample 2
# unidimensionality ND
fa3_model_scaleND <- "
fa3_scaleND =~ FEQ_04 + FEQ_05 + FEQ_07 + FEQ_09 + FEQ_11 + FEQ_12 + FEQ_24 + FEQ_27 + 
              FEQ_36 + FEQ_37
"

s2.fa3.cfa_scaleND <- cfa(fa3_model_scaleND, data = feq_sample_2, estimator = "MLM", mimic = "Mplus")

summary(s2.fa3.cfa_scaleND, fit.measures = T, standardized = T)


# reliability ND
s2.fa3_scaleND <- feq_sample_2 %>% select(FEQ_04, FEQ_05, FEQ_07, FEQ_09, FEQ_11, FEQ_12, 
                                          FEQ_24, FEQ_27, FEQ_36, FEQ_37)

set.seed(42)
ci.reliability(data = s2.fa3_scaleND, type="hierarchical", conf.level = 0.95,
               interval.type="bca", B=1000)
alpha.fa3_scaleND = alpha(s2.fa3_scaleND)
print(alpha.fa3_scaleND, digits = 3)